Comparing Covid Tests and Covid Positives using mapping data

Through mapping both covid tests and positives adjusted for population, we hope to fine if any relation between the two exists. It may be possible that more tests mean more covid.

Load packages

The following packages are used:

knitr– load datasets

tidverse– mutate data

janitor– clean names

lubridate – work with dates

here– give file locations

library(knitr)
library(tidyverse)
library(janitor)
library(lubridate)
library(here)

Additional packages being used

tmap– for mapping purposes

sf– joining spacial data and generating plots

tidycensus– loading census data

library(tmap)
library(sf)
library(tidycensus)

Load testing data

The test data is loaded into R, and mutated to become non-cumulative

df1=readxl::read_excel(here("submissions","bike_share","data","tot_tests_neighborhood.xlsx")) %>% 
  clean_names() %>%
  #mutate(date_reported= as_datetime(date_reported))
  separate(date_reported,into=c("date","hr"),sep = " ")%>%
  mutate(date_report=as_date(date))%>%
  group_by(neighborhood)%>%
  mutate(across(total_tests, ~ .-c(0,lag(.)[-1])))

Load Mapping Data

Loads in the Health neighbourhood geojson file.

neigh=st_read(here("submissions","bike_share","data","dc_neigh.geojson")) %>% clean_names()

Reading layer `DC_Health_Planning_Neighborhoods' from data source 
  `D:\KALTENBACH_ENGINEERING\Master\CU\R_Workspace\Portfolio_DS\submissions\bike_share\data\dc_neigh.geojson' 
  using driver `GeoJSON'
Simple feature collection with 51 features and 8 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -77.11976 ymin: 38.79165 xmax: -76.9094 ymax: 38.99556
Geodetic CRS:  WGS 84

class(neigh)

[1] "sf"         "data.frame"

Select and clean data

Cleans testing data so it can be joined with the health neighborhood spacial data

df_tests=df1 %>%
  filter(as_date(date_report) == "2021/10/27") %>% 
  separate(neighborhood,into=c("code","name"),sep = ":") %>% 
  mutate(code=case_when(
    code=="N35" ~"N0",
    TRUE ~ code
  ))

Join mapping data with test data

Joins data and generates plot

neigh2=left_join(neigh,df_tests,by=c("code")) 
tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(neigh2) +tm_polygons("total_tests",alpha=.5)

Load positive cases data

Reads in positive coivd cases data and cleans it so it can be joined to the above datasets

df_c=readxl::read_excel(here("submissions","bike_share","data","neigh_cases.xlsx"),
                            col_types = c("date", "text", "numeric")) %>% 
  clean_names() 
df_cases=df_c %>%
  filter(as_date(date) == "2021/10/27") %>% 
  separate(neighborhood,into=c("code","name"),sep = ":") %>%
  mutate(code=case_when(
    code=="N35" ~"N0",
    TRUE ~ code
  ))

Join covid data

neigh3=left_join(neigh2,df_cases,by=c("code")) 
tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(neigh3) +tm_polygons("total_positives",alpha=.5)

Graph data

 tm_shape(neigh3)+tm_polygons(c("total_positives","total_tests"),alpha=.5)

Get census Data

Census data is used to normalize both tests and COVID data for population

census_api_key("2c8b9d5c4902b7efb4e1f98b2c23692cb1b73e95")

To install your API key for use in future sessions, run this function with `install = TRUE`.

v20 = load_variables(2018,"acs5")
df_cencus=get_acs(geography = "tract",
                  variables=c("pop"="B01001_001"),
                  state="DC",geometry=TRUE,year=2018) 

Getting data from the 2014-2018 5-year ACS
Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
Using FIPS code '11' for state 'DC'

df_cens=df_cencus %>% select(-moe) %>% spread(variable,estimate) 
df_cens_adj=df_cens %>% st_transform(4326)

Join Census data to mapping data

df_j=st_join(df_cens_adj,neigh3,largest=TRUE)

Warning: attribute variables are assumed to be spatially constant throughout all geometries

Summarize population for each neighbourhood

df3=df_j %>% select(pop,code) %>%
  group_by(code) %>%
  summarise(pop_n=sum(pop))
            
            

Join population data and generate results

Removed any extreme outliers and generated map to compare testing rates and covid rates

df4=left_join(neigh3,df3 %>% st_set_geometry(NULL))

Joining, by = "code"

df4=df4 %>% mutate( covid_rate=total_positives/pop_n, test_rate=total_tests/pop_n)
df4 %>% filter(!(code %in% c("N0","N15","N24"))) %>% tm_shape()+tm_polygons(c("covid_rate","test_rate"))